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We numerically investigate the thermally unstable accretion disks around 
black holes. We adopt an evolutionary viscous stress equation to replace the 
standard alpha-prescription based on the results of two MHD simulations. We 
hnd a kind of interesting oscillations on some running models in limit-cycle out¬ 
burst state. The oscillations arise near the inner boundary and propagate radi¬ 
ally outwards. We deem that they are the trapped p-mode oscillations excited 
by sonic-point instability. We directly integrate the local radiation cooling fluxes 
to construct the mimic bolometric light-curve. We hnd a series of overtones be¬ 
side the fundamental harmonic on the power spectra of mimic light-curves. The 
frequency of the fundamental harmonic is very close to the maximum epicyclic 
frequency of the disk and the frequency ratio of the fundamental harmonic and 
overtones is a regular integer series. We suggest that the code for ray-tracing 
calculation must be time-dependent in virtual observation and point out the ro¬ 
bustness of the black hole spin measurement with high frequency QPOs. 

Subject headings: accretion, accretion disks — black hole physics — gravitation 
— relativistic processes — X-rays: bursts 


Introduction 


High-frequency quasi-periodic oscillations (HFQPOs) have been observed in some black 
hole (BH) X-ray binaries, which only appear in "steep power law" state in high luminosity 
(L > O.lLEdd) and are in range of 40 to 450Hz. These frequencies are comparable to the 
orbital frequency of innermost stable circular orbit (ISCO) of a stellar-mass BH. It is believed 
that the mechanism behind them is closely related to the dynamics of inner regions of BH 
accretion disks (see Remillard fc McClintock 2006 : Kato et al 2008 : Belloni et al 2012 . for 


reviews 


The pioneering work of Katol 1 1978 ) studied the effect of viscosity on the disk stability 
and found the sonic-point instability (its original name is pulsational instability), which is the 
excitation process of so called p-mode oscillations in transonic accretion flows (this oscillation 
has the other name, inertial-acoustic oscillation, but we only call it p-mode oscillation in 
this paper). The mechanism of this instability is anal ogous to the e-mechanism in stellar 
pulsations. If the viscous parameter a of a-prescription fjShakura fc Snnvaev 1119731) is larger 
than a critical value, the instability will arise from the phase relation between the viscous 
heat generation and oscillation, i.e. the viscous heat generation increase (de crease) in the 
compressed (expanded) phase of disk oscillation (jKatolll978l: iKato et allll988f) . 
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A series of works flMatsumoto et allllQSSl. Il989[ IChen &: Taainlll995l: iMiranda et a]|l2014i) 
were dedicated to study the role of viscosity on HFQPOs. They all found the oscilla- 
tion, whose fr e quenc y is close to the maximum of epicyclic frequency ^max- Additionally, 
Miranda et all (120141) also found many overtones, whose frequencies are close to the integer 
multiples of Kmax- 

Though the famous a-prescription has been used extensively since 1973, it may be too 
simple to closely accord wi th actual ac c retion flows. The magnetohydrodynamic (MHD) 
shearing box simulations of lHirose et all (1200911 implied that there is certain time-delay be¬ 
tween the viscous stress and tot al pressure instead of the instantaneity introduced by a- 
prescription. The other work of iPenna et all (120131) . which is based on a few relativistic 
MHD global simulations, pointed out that the parameter a is a function of radius but not 
constant in the inner disk region. More than one decade ago, lYamasaki &: Kato (jl996[) 


analyzed effects of viscous time-delay on the p-mode oscillations and they found that it has 
a little of inhibition on the growth of oscillations but excitation still exists for high viscosity 


cases. 


Our work presented in this letter is also dedicated to study the role of viscosity on QPOs. 
The distinction of our model is the evolutionary viscous stress equation (instead of the a- 
prescription), which is constructed to mim ic ef fects of viscous tiin e-delay and inconstant a 
basing on the results of lHirose et all (120091) and iPenna et all (120131) . 


2. Equations 


In this letter, we consider the axisymmetric relativistic accretion flows around black 
holes. The Boyer-Lindquist coordinates t, r, 9, cj) are used to describe the space-time. The 
governing equations of accretion flows are written as following. 
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Where S, 4^, £, H, U, T, and Sr,p are the surface density, radial velocity (measured in 
the corotating frame), angular momentum per unit mass (£ = U(f,), half thickness of disk, 
vertical velocity of the surface, local temperature of accreted gas, and viscous stress (only 
r^-component is non-vanishing), respectively. These equations are all derived from the 
conservat i on of t he stress-energy tensor, and almost the same as those in our previous paper 
Xne et ai (2011) (see that paper for the detailed derivations as well as definitions of 7 , A, 
A, and etc). 

As mentioned in previous section, we adopt an additional evolutionary stress equation 
to describe the viscosity instead of the a-prescription. This equation can be written as 

.dSrcj, 
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where the factor nr* is the practical viscous time-delay, which is scaled with the typical 
delay r* by parameter n; S*^ is the expected stress by turbulence. The definitions of r* and 


S*^ can be written as 
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If a = const a nd n —?• 0, equat ion ()7|) will reduce to = S*^, which is the same as the a- 
prescription in IXue et all (120 111) . It means that equation (|7j) contains the a-prescription as a 
trivial case. Equation (IT^ determines the dependence of a on rad ius r, BH mass M and spin 
a. The radial factor, including the exponent 6 , was suggested bv iPenna et all (120131) . Under 
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this profile, a is almost a constant ckq in enter disk region with large r and increases to higher 
valne radially inwards. We set ao = 0.1 and fix the mass supplying rate Mout = O.OOMEdd 
(see table [T] for the definition of Eddington accretion rate, MEdd) at the outer boundary in 
this work. In practice, these settings are sufficient to make the disk thermally unstable and 
we indeed observe the limit-cycle outburst from running code. 

After dehning S'^^, the viscous heating rate F~^ in equation ([6]) is redehned as 


^yll/2AV2 OQ 
r — -/-^ 

CJ'Y' 


which would be reduced to the one in IXue et all (120111 ) when n —)■ 0. 


(13) 


3. Instabilities 


We update our previous code established in IXue et all (120llh an d ru n it for eleven 
numer ical models, whose parameters are listed in tabled] iLin et all (120111) and ICiesielski et al 
(120121) studied the impact of viscous time-delay on the thermal instability and found that it 
is not remarkable for small enough time-delay. Indeed, we observe the expected limit-cycle 
outbursts on all models with the time delay parameters n in the range 0 to 4. Therefore, 
the time-delay implemented in our code does not affect the limit-cycle outbursts. 

Among these models, S8 and S30 are two typical ones. S8 is a disk around a non¬ 
spinning black hole as well as S30 around a fast-spinning one. In figure dl we show the 
bolonietric light curves for these two models respectively. For each point on light-curves, the 
disk luminosity is made by integrating the local radiation cooling fluxes, and the sampling 
time corresponds to the frequency lO^Hz, which is enough to reveal any harmonic with 
frequency lower than 5 x lO^Hz in the power spectral density (PSD). Due to the difficulties 
in hydrodynamical calculation, we only obtain one outburst light curve for S8 and 0.7s-long 
luminosity ascending light curve for S30. However, they are long enough for the calculation 
of PSDs. In figure [21 we show the relevant PSDs for S8 and S30. The fundamental frequency 
(the lowest frequency of harmonics) is ~ 74.9Hz for S8 and ~ 285.6Hz for S30, which are 
both close to 71.3Hz and 300Hz, the respective maximal epicyclic frequencies, which are 


Kato 

to 

o 

o 

1 -^ 

Kato et al 

2008 


and a detailed relativistic analysis is ongoing by 
our colleague Jiff Hor ak). The spect r um of axisymmetric (m = 0), horizontal p-modes was 
recently computed by iGinssani et all (1201411 , who show that in addition to a discrete set of 
lower frequency modes which are trapped in the inner disk, there are modes of frequency 
very close to the maximal epicyclic frequency in which the oscillation is transmitted to the 
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outer disk. In table [H we also list the fundamental frequencies observed from the other 
oscillating models, which are all close to the theoretical values though they vary in a narrow 
frequency range. 


In hgure [3l we show an oscillating cycle of model S8 for V (upper two panels) and Sr^ 
(lower two panels). The oscillations arise near the inner boundary and propagate outwards. 
In the hgure, one can follow the motion of individual wavelets. The negative V denotes 
inhow and negative gradient of V corresponds to the compression (inhow speed of inner is 
slower than outer) as well as positive gradient to the expansion (inhow speed of inner is 
faster than outer). Thus, any wavelet on the l/-prohles can be divided into the compressed 
wave-front and expanded wave-rear regions. For example, the left-most wavelet on initial V- 
prohle (red curve in the upper-left panel) can be divided into the wave-front (between dashed 
and dash-dotted lines) and wave-rear (between solid and dashed lines) regions. The relevant 
variation of which is proportional to the viscous heating, is showed in the lower-left 
panel. Sr^ (as well as the viscous heating) monotonously increases in compressed wave-front 
(see the red curve between dashed and dash-dotted lines) but it unceasingly increases after 
the maximal compression (takes place at the dashed line) and then decreases in the expanded 
wave-rear. This non-monotonicity of in wave-rear is due to the time-delay contained by 
equation [71 We deem that this is the phase relation between the viscous heat generation 
and disk os c illation, which im plies the arise of sonic-point instability for p-mode oscillation 
fjKatolll978l: iKato et allll988[) . Thus, we also deem that the oscillations actually arise near 
the sonic-point, which is included in our computational domain and very close to the inner 
boundary. 


Results and Discussions 


In table [H there are the other nine models with different viscous settings but with the 
same BH mass, BH spin and mass supply rate (fixed accretion rate at the outer boundary) 
as S8. We only observe oscillations on the models with non-vanishing delay and large enough 
constant a (> 0.3) or a-prohle. In fact, the effective value on the a-pro£le increases inwards 
from constant 0.1 in outer disk region to the maximum 0.45 at the inner boundary. Thus, 
the impact of a-prohle is similar to the large constant a while it just becomes large enough 
near the origin of instability, sonic-point. The facts of these oscillating models imply that 
the large a (at least near sonic-point) is a necessary condition for the arise of sonic-point 
i nstab ility as well as p-mode o s cillat ion, which is consistent with the analysis of iKato et al 
(1988) and Yamasaki &: Kato ( 1996 ). 


On the other view, we note the effect of delay on the appearance of oscillations. The 
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oscillating models with n = 1 (S7 and S8) have the oscillations during the whole outburst. 
The other oscillating models with n ^ 1 (SI, S2 and S17) lose the oscillations in different 
luminosity stages. The models with n = 0 (S3, S4, and S16) have no any oscillation though 
they have large enough a, while the disappearance of the oscillations on models S5 and S15 
is due to the small a. These facts imply that r* may be a favorable delay for oscillation 
excitation on the disk around a IOMq Schwarzschild BH under the large a. The last model 
in table [H S30 is a special case for a fast-spinning Kerr BH, which has another favorable 
delay 4r*. This may imply the dependence between the viscous time-delay and BH spin. 


Focusing on the luminosity of the oscillating models, we observe oscillations only in the 
limit-cycle outburst state (L > 0.2LEdd) when the inner disk region has switched to slim disk 
mode. On the contrary, there is no any oscillation observed in the limit-cycle quiet state 
(L ~ O.OlLEdd)- This is consistent with the HFQPO observations, but cannot be compared 
with the sonic-point instability theory w hich does not discri minate between accretion rates. 
Recently, the shearing box simulation of lHirose et all (1201411 implied that the effective a is 
en hanced by th e verti cal convection during the outburst, which is similar to the conception 
of iMilsom et all (1199411 . Thus, larger a required by the p-mode oscillation may be caused by 
the outburst, explaining why HFQPOs are observed only in high luminosity state. 


Beside the fundamental harmonic, there are many overtones in both of the two spectra in 
figure 121 The frequency ratio of fu ndamental harnionic and its overtones is a regular integer 
series, which is also observed by Miranda et all (1201411 in their 2D-simulations (in radial 
and azimuthal dimensions) when the a xisymmetric p-mode become dominant for large a. 
However, no overtones are observed bv I Chen fc TaamI (1199511 . Perhaps, this is due to their 
adopt ing of the viscous prescription with Sr(f, oc Pgas instead of Sr,p oc ptotai in Miranda et al 
(1201411 and our models. These interesting overtones may be potentially useful for explaining 
the observational QPO pairs, which are always in a specific integer frequency ratio. 


Subsequently, the virtual observation from our numerical results will be logical and 
interesting. However, it would require more careful treatment. For example, the effective 
time-delay of arrival caused by the gravitational bending on the trajectories of emitted pho¬ 
tons and the large observational view-angle, and effective blocking caused by the gravitational 
red-shift and the other shields. 


As a rough evaluation for the effective time-delay of arrival, one can consider the obser¬ 
vation from an almost edge-on disk, on which the photons emitted from two locations apart 
from distant Ar = 203M at the same observer’s time (t) will arrive at the observer with the 
rough delay At = 0.01s comparable with the period of 71.3Hz (our results presented in fig¬ 
ures [T] and I2] can be roughly regarded as from the face-on disks). The distant Ar = 203M is 
comparable with the radius of outward wavelets propagating area, so it is possible to change 
















observer’s final view on the oscillation power spectrnm. This also implies that the code for 
ray-tracing calculation on virtual observation must be time-dependent. 

In order to roughly demonstrate the effective blocking, we calculate the light-curves 
without the radiation contribution from different inner cutting regions for the same model 
S8 and we show the relevant power spectra in figure 01 It is remarkable that the fundamental 
harmonic (inside the rectangle in all four panels) cannot be easily removed from PSDs 
because of the outward propagation of the oscillation from sonic-point. It implies that the 
measurement of BH spin with HFQPO will be very robust even in a case when modulation 
of the innermost disk is not visible. 


So far, we have found some features of our model fortunately coinciding with the coun¬ 
terparts of HFQPO observations. All of these are only due to the adopting of a special 
evolutionary stress equation in our model, which mimics the viscous features induced from 
MHD simulations. However, our model still lacks some abilities for capturing various com¬ 
plicated features associated with the real accretion flows around BHs. In fact, it is almost 
impossible for seeing the oscillations during the whole outburst state (like those in figure 
[T]) in real accretion flows. It is because the effective a and viscous time-delay, determining 
the appearance of oscillations, is turbulent stochastic in real accretion flows though their 
expected values can be determined by certain laws. Further more MHD simulations on the 
mean behaviors of effective turbulent viscosity is necessa r y for i mproving our und erstanding 


on accretion process around BHs (e.g. iHirose et alll2009l. I2014J: iPenna et all 120131) . It is also 


im possible for cap t uring the resonance between the radial and vertical oscillations suggested 
bv iKluzniak et all (120041) . because our model is a vertically integrated model and there is 
only radial dependence reserved. The further analytic works and MHD simulations on the 
roles of radial and vertical oscillations are necessary for explaining the observational HFQPO 
pairs though our model has the intrinsic multi-frequency feature. In fact, the thermal insta¬ 
bility required by our model is the theoretical one whose deviations from the observations 
have been found a decade ago flOierlinhski &: Donel 120041) . While the existence of thermal 
instability on the BH accretion disks is still an open issue at present, we only can adopt 
this theoretical thermal instability to produce the high luminosity outburst required by the 
oscillations in our model. We believe the plausibility of oscillations observed on our models 
though it is "dancing" on a poor-quality stage. 


This work was supported by the National Natural Science Foundation of China un¬ 
der grants 11233006 and 11373002, Polish NCN grant UMO-2011/01/B/ST9/05439 and 
2013/08/A/ST9/00795, and Czech ASCRM100031242 CZ.1.07/2.3.00/20.0071 Synergy (Opava) 
project. 
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Table 1. Model Sequences 


ID 

n 

OL 

M/Mq 

ajM 

M/MEdd^ 

Fundamental 
Frequency [Hz] 

Max Epicyclic 
Frequency [Hz] 

S8 

1 

Eq. |(T2)| 

10 

0 

0.06 

71.64-76.732 

71.3 

S15 

1 

0.1 

10 

0 

0.06 

No Osci. 

71.3 

S5 

1 

0.15 

10 

0 

0.06 

No Osci. 

71.3 

S7 

1 

0.3 

10 

0 

0.06 

62.37-68.7 

71.3 

S4 

0 

0.3 

10 

0 

0.06 

No Osci. 

71.3 

S3 

0 

0.45 

10 

0 

0.06 

No Osci. 

71.3 

S16 

0 

Eq. |[T2t 

10 

0 

0.06 

No Osci. 

71.3 

S2 

0.5 

Eq. Hill 

10 

0 

0.06 

74.77-80.64 (AD)^ 

71.3 

S17 

0.75 

Eq. (HU 

10 

0 

0.06 

74.27-78.5 (AD) 

71.3 

SI 

2 

Eq. (HU 

10 

0 

0.06 

69.77-71 (P) 

71.3 

S30 

4 

Eq. (HU 

7.02 

0.947 

0.06 

285.64 4 

300 


^Eddington accretion rate Msad = ~ 2.23 x 10®-j^-(gs ^). 

^The variation range of observed fundamental frequencies is given. 

^The characters inside the parentheses denote the luminosity stages in which the oscillations 
disappear. For example, (AD) denotes the oscillations only appear in the luminosity plateau 
stage (disappear in ascending and descending stages). 

^We only have the data for the luminosity ascending stage. Due to the lack of data, we cannot 
observe any remarkable variance on its fundamental frequency. 
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Time [second] Time [second] 


Fig. 1.— Light-curves of models S8 (left) and S30 (right). There are a few subplots to reveal 
detailed views of surrounding arrow points. 



Fig. 2.— PSDs of models S8 (left) and S30 (right). 
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Fig. 3.— One oscillation cycle of the radial velocity and viscous stress from model S8. The 
different color lines denote a serial snapshotting times, which are all scaled with the period 
of this cycle. The radius r is scaled with BH mass M because we take G = c = 1. Thus the 
Schwarzschild radius is 2M. 
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Fig. 4.— PSDs of different blocking cases. From upper to lower, the radii of cutting regions 
Tcut increase for the same model S8. 

























